* Reset settings and initialize log file
launch, path("share/couples")

*-------------------------------------------------------------------------------
* Price and Wasserman (2024), "The Summer Drop in Female Employment"
*
* Description: Show joint employment rates among heterosexual married couples.
*-------------------------------------------------------------------------------


* Prepare data and run models
*-------------------------------------------------------------------------------

if "$estimate" != "0" {
	* Load data
	gzuse "$basepath/data/derived/cps_bms_sample.dta.gz", clear
	keep pid hid tm month tmspline* weeks wtfinl female spouse_present youngest pernum sploc emp

	* Restrict to households with (1) exactly one adult man and one adult woman
	bysort hid tm (female): keep if _N == 2 & female[1] != female[2]

	* Restrict to households where the two adults are married to each other and both are present
	bysort hid tm (female): keep if spouse_present[1] == 1 & spouse_present[2] == 1
	bysort hid tm (female): assert (pernum[1] == sploc[2]) == (pernum[2] == sploc[1])
	bysort hid tm (female): keep if pernum[1] == sploc[2]
	drop pid spouse_present pernum sploc

	* Verify that the age of the youngest child is constant within couples
	bysort hid tm (female): assert youngest[1] == youngest[2]

	* Weight households by average of the spouses' sampling weights
	bysort hid tm: replace wtfinl = 0.5 * (wtfinl[1] + wtfinl[2])
	bysort hid tm: replace wtfinl = wtfinl[1]

	* Measure combinations of male and female employment
	reshape wide emp, i(hid tm) j(female)

	* Create indicators for joint employment status
	gen joint_ee = (emp0 == 1 & emp1 == 1)
	gen joint_en = (emp0 == 1 & emp1 == 0)
	gen joint_ne = (emp0 == 0 & emp1 == 1)
	gen joint_nn = (emp0 == 0 & emp1 == 0)

	* Create a set of observations for parents with young school-age children
	expand = 2 if youngest == 2
	bysort hid tm: gen byte group = _n

	* Compute aggregate employment rates
	gcollapse (mean) joint_* (rawsum) wtfinl (first) month tmspline* weeks [pw = wtfinl], by(group tm)
	tsset group tm

	label define group_lbl 1 "All married (and co-residing) couples", replace
	label define group_lbl 2 "Youngest child aged 6-12", add
	label values group group_lbl

	* Run regressions
	foreach g of numlist 1 2 {
		foreach j in "ee" "en" "ne" "nn" {
			* Run specification
			quietly ivreg2 joint_`j' ib5.month tmspline* weeks if group == `g' [aw = wtfinl], bw($bandwidth) robust small
			process_estimates, path("couples") model("g`g'_`j'")
		}
	}
}


* Prepare estimates
*-------------------------------------------------------------------------------

* Load estimates into memory
load_estimates, path("couples")

* Prepare coefficient labels for each month
make_coeflabels


* Plot seasonality in joint employment and presence at work
*-------------------------------------------------------------------------------

* Prepare background shading
clear
set obs 14
gen rlow = -4
gen rupp = 4
gen rval = _n - 0.5

* Plot estimates
foreach g of numlist 1 2 {
	if `g' == 1 local glbl "All married (and co-residing) couples"
	if `g' == 2 local glbl "Youngest child aged 6-12"

	#delimit ;
	coefplot
		(g`g'_ee, offset(-.15) recast(connected))
		(g`g'_en, offset(-.05) recast(connected))
		(g`g'_ne, offset(+.05) recast(connected))
		(g`g'_nn, offset(+.15) recast(connected)),
		coeflabels(`coeflabels')
		title("{bf:`glbl'}")
		xtitle("")
		xlabel(, alternate)
		ytitle("")
		yscale(range(-4 4))
		ylabel(-4(1)4)
		addplot(
			rarea rupp rlow rval if inrange(rval, 0.5, 5.5), color($ltgs) plotregion(margin(t=0 b=0)) below ||
			scatteri 0 0.5 0 13.5, recast(line) color(black) below)
		vertical
		legend(rows(2) size(*0.8) bmargin(0 0 0 2) keygap(1.5) order(
			4 "Both employed"
			6 "Only husband employed"
			8 "Only wife employed"
			10 "Neither employed"))
		rescale(100);
	#delimit cr

	tempfile g`g'
	graph save "`g`g''"
}

* Create pdf file
grc1leg "`g1'" "`g2'", l1title("      Difference relative to May (p.p.)", size(*.8)) imargin(vsmall)
nicepdf "$basepath/output/couples.pdf", panels(2) indirect replace

* Report point estimates
estimates replay g1_ee

* Close the log file
unlaunch
